ATTITUDE ESTIMATION OR QUATERNION ESTIMATION? 

F. Landis Markley* 

The attitude of spacecraft is represented by a 3x3 orthogonal matrix with unity 
determinant, which belongs to the three-dimensional special orthogonal group 
SO(3). The fact that all three-parameter representations of SO(3) are singular or 
discontinuous for certain attitudes has led to the use of higher-dimensional 
nonsingular parameterizations, especially the four-component quaternion. In 
attitude estimation, we are faced with the alternatives of using an attitude 
representation that is either singular or redundant. Estimation procedures fall into 
three broad classes. The first estimates a three-dimensional representation of 
attitude deviations from a reference attitude parameterized by a higher- 
dimensional nonsingular parameterization. The deviations from the reference are 
assumed to be small enough to avoid any singularity or discontinuity of the 
three-dimensional parameterization. The second class, which estimates a higher- 
dimensional representation subject to enough constraints to leave only three 
degrees of freedom, is difficult to formulate and apply consistently. The third 
class estimates a representation of SO(3) with more than three dimensions, 
treating the parameters as independent. We refer to the most common member of 
this class as quaternion estimation, to contrast it with attitude estimation. We 
analyze the first and third of these approaches in the context of an extended 
Kalman filter with simplified kinematics and measurement models. 

INTRODUCTION 

Real-time spacecraft attitude estimation generally employs an Extended Kalman Filter (EKF) [1, 2]. 
Although the 3x3 orthogonal attitude matrix is the fundamental representation of the spacecraft’s attitude, 
the orthogonality requirement imposes six constraints on its nine elements, reflecting the fact that the 
special orthogonal group SO(3) of rotation matrices has dimension three. Therefore, employing the nine 
elements of the attitude matrix as components of a state vector in a Kalman filter leads to some complexity 
in enforcing the constraints [3, 4], Most EKFs use lower-dimensional parameterizations of SO(3), and 
several have used minimal three-dimensional parameterizations [5-7]. However, it is a well-known fact that 
all three-parameter representations of SO(3) are singular or discontinuous for certain attitudes [8]. This has 
led to the use of higher-dimensional nonsingular parameterizations in EKFs, especially the four-component 
quaternion [9-13]. Reference 12 presents an overview of Kalman filtering for spacecraft attitude estimation, 
emphasizing the quaternion representation, with a complete list of references through 1981. 

The quaternion has the lowest dimensionality possible for a globally non-singular representation of SO(3), 
but it still has one superfluous degree of freedom. Thus we are faced with the alternatives of using an 
attitude representation that is either singular or redundant. Various strategies to avoid or evade this dilemma 
have been proposed and analyzed [14-18], and we will consider two in this paper. Our preferred strategy, 
which we refer to as the Multiplicative EKF (MEKF), uses a nonsingular representation for a reference 
attitude and a three-component representation for the deviations from this reference. This method, which was 
implicit in some early EKFs for estimating spacecraft attitude [9-1 1], has been discussed in detail in 
Reference 18. The second strategy we discuss treats the four components of the quaternion as independent. 
We also mention, but we do not analyze in detail, methods that estimate a constrained four-component 
quaternion with three degrees of freedom. These methods have been treated extensively elsewhere [12-17]. 
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We begin with a review of the attitude parameterizations we will use [19, 20], The MEKF is then discussed 
in a model with simplified kinematics and measurement processing, in order to bare its essentials. We next 
turn to the concept of a quaternion probability distribution function, which is not as straightforward as one 
might expect. This is followed by an analysis of three versions of an EKF for a quaternion with four 
independent components, using the same models as for the MEKF, and by a summary of our conclusions, 

QUATERNIONS 

A quaternion is a four-component object with a three-vector part and a scalar part 


q 


The four components of q , known as the Euler symmetric parameters or the Euler-Rodrigues parameters, 
first appeared in a paper by Euler [21] and in unpublished notes by Gauss [22]; but Rodrigues’ classic paper 
of 1840 first demonstrated their general usefulness [23]. Hamilton introduced the quaternion as an abstract 
mathematical object in 1844 [24], but there is some question as to whether he correctly understood its 
relation to rotations [25], The attitude matrix is generally written as a homogeneous quadratic function of 
the components of a quaternion, 

A( 5) = (<74 -N 2 ) / 3x3 + 2qq T -2q 4 [qx], (2) 

where / 3x3 is the 3x3 identity matrix, and the cross product matrix is 


[q x ] s 


0 

<7 3 

L-<7 2 


<73 <72 

0 -<?, 

<7i 0 


(3) 


We see from Eq. (2) that q and -q represent the same attitude, and that this equation gives an orthogonal 
matrix only if the quaternion q has unit length; 


In this case we can write 


fef-tf+rf-i. 


<7 


esin(0/2) 

cos(0/2) 


(4) 

(5) 


where e and <p are the Euler axis and angle of rotation, respectively. We follow References 12 and 20 in 
writing the quaternion product as 


P 


_ m+^p-pxq 

<7 S 

R4<?4 -p-q 


( 6 ) 


This differs from the historical multiplication convention [19, 24], denoted by pq without an infix 
operator, by the sign of the cross product in the vector part. The two products are related by p®q=qp. 
The convention adopted here has the useful property that 

A(p)A(q) = A{p®q). (7) 

With the historical convention, the quaternion ordering on the right side of the above equation would be the 
reverse of the order on the left side. With either convention, though, the product of two quaternions is 
bilinear in the elements of the component quaternions, a property shared with the direction cosine matrix 
but no other attitude representation. In fact, Eq. (7) means that the rotation group and the quaternion group 
are almost isomorphic, the qualifier “almost” owing to the 2: 1 nature of the mapping [26], 
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The quaternion obeys the kinematic equation 


< 7 = 2 




(8) 


where co is the angular velocity vector in the body frame, and we use a boldface symbol with an overbar to 
denote a quaternion with zero scalar part. 


MULTIPLICATIVE EXTENDED KALMAN FILTER 


The MEKF represents the attitude as the quaternion product 


q=8q(a)®q ref , (9) 

where q ref is some unit reference quaternion and 8q{ a) is a unit quaternion representing the rotation from 
q rrf to the true attitude q , parameterized by a three-component vector a. Although several choices for a are 
possible [18], in this paper we choose it to be two times the vector of Rodrigues parameters [19, 20, 23], 


a = 2<5q/(5<7 4 = 2etan(0/2) . 


( 10 ) 


The factor of two is chosen so that |a| = when the rotation angle between q and q nf is small. The 


inverse of Eq. (10) is 


Sq(a) = 


1 

M 

y/ 4 + W 2 

_2_ 


( 11 ) 


The three Rodrigues parameters are often referred to as the Gibbs vector, since they were arrayed by Gibbs in 
his “vector semitangent of version” [27]. They are infinite for 180° rotations, which is undesirable for a 
global representation of rotations; but they provide an excellent representation of rotations less than 180°. 
Since Eq. (10) maps q and -q to the same point in three-dimensional Euclidean space E3, these parameters 
provide a 1:1 representation of the rotations less than 180° onto E3, which provides support for Gaussian 
distributions with infinitely long tails. 


The two attitude representations Sq and q ref in Eq. (9) are clearly redundant. The basic idea of the MEKF is 
to compute an unconstrained estimate of the three-component a while using the correctly normalized four- 
component q ref to provide a globally nonsingular attitude representation. Given an estimate a, where a 
caret will denote the expectation of a random variable, Eq. (9) says that the corresponding estimate of the 
true attitude quaternion is 8q( a) ® q re/ . We remove the redundancy in the attitude representation by 
choosing the reference quaternion q re/ so that a is identically zero. Since <5*7(0) is the identity quaternion, 
this means that the reference quaternion is the best estimate of the true quaternion. Note that the reference 
quaternion in the MEKF is not considered to be a random variable, and its estimate is not an expectation. 
The fundamental conceptual advantage of the MEKF is that q ref is a unit quaternion by definition. 


The identification of q ref as the attitude estimate means in turn that a is a three-component representation 
of the attitude error, the difference between the truth and our estimate. This provides a consistent treatment 
of the attitude error statistics, with the covariance of the attitude error angles in the body frame (in radians 
squared) represented by the covariance of a. 


An alternative formulation, which has some advantages, reverses the order of multiplication in Eq. (9) so 
that a represents the attitude errors in the inertial reference frame rather than in the body frame [28]. It is 
also possible to represent the reference attitude by a reference attitude matrix A ref rather than by a reference 
quaternion [9, 28], This requires more parameters, but may save computations if the attitude matrix is 
explicitly required. An argument in favor of the quaternion is that it is easy to restore normalization that 
may be lost due to numerical errors, while restoring the orthogonality of A ref is nontrivial. Gray has argued 
that this argument is not compelling if reasonable computational care is taken [28]. 
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Kinematics 


Since q ref is a unit quaternion, it must obey a kinematic equation of the form 

( 12 ) 

where <& ref has the obvious interpretation as the angular velocity of the reference attitude. We now show 
how ( 0 rf/ is determined by the requirement that a be identically zero, which is the condition that the 
reference attitude be the optimal attitude estimate. Computing the time derivative of Eq. (9), using Eqs. (8) 
and (12), gives 

^®q=S7j®q ref + ^8q®(5 ref ®q ref . (13) 

Substitute Eq. (9) for q on the left side, right-multiply the entire equation by the inverse of q ref , and 
rearrange to get 

8q = ^(<5®8q-8q®<5 ref ). (14) 

Note that co rc/ is not a random variable. Substituting Eq. (14) into the time derivative of Eq. (10) gives 


a = (/jxj + 1 aa 7 ) (to -oo ref ) - \ (o + co rc/ ) x a, 

after some straightforward quaternion algebra. We now assume the simple kinematic model 

a(t)=&( t) + n m (t), 

where to (t) is the nominal angular velocity and (t) is a white noise process, satisfying 

E{ n a (t)nl (t ‘ ') } = Q S(t - f ') = cl J 3x3 S(t-t'), 
with 8{t - 1') denoting the Dirac delta function. Linearizing Eq. (15) in a and gives 

a =(0 -03 ref + Fa + Gn B 

with 

Fs -![(© + © r(/ )x] 
and 


(15) 

(16) 

(17) 

(18) 

(19) 

( 20 ) 


Measurements 

We assume discrete measurements with covariance R. For ease of analysis, though, we will assume that the 
measurement interval r is so much shorter than any filter time constant that we can treat the measurements 
as continuous with covariance Rx [1], The attitude error estimate a conditioned on the measurements is 
then given by the expectation of Eq. (18) augmented by a measurement term [1], 

a = ® ~ ® ref + Fa + P a H T M (tR)-' [ z - h(a)] , (21) 


where the subscript M denotes the MEKF and the 3x3 covariance matrix of the attitude error vector satisfies 

Pa = PPa + P«P T +GQG t - P a H T M (T/?r H m P a . (22) 

We model the measurement of a vector \ B in the body coordinate frame by 

z ~ h w (a) + a, = v g -A(q) + n z , (23) 
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where v, is the corresponding vector in an inertia! frame and n. is Gaussian measurement noise. To first 


order in the attitude error vector a, we have 

A(q) = A{Sq®q rtf ) = A(Sq)A(q rrf ) = (/ 3x3 -[a x])A(q rtf ), (24) 

where the last step follows from Eqs. (2) and (1 1). Substituting into Eq. (23) gives 

M a ) = v s -(^* 3 -[ ax ]M(g re/ )v, =v B -v M - v M x a (25) 

with v M denoting the measurement in the body frame predicted by the reference quaternion: 

= A(q rcf )\,. (26) 

It follows that the measurement sensitivity matrix is 

h m s 3h w /aaj. =-[v w x] . (27) 

If we assume that the vector measurement errors are isotropic, as is appropriate for measurement of a full 
vector such as the magnetic field, we have 

R=c 2 Jw (28) 

Making appropriate substitutions into Eq. (21) gives 

a = ® - <*»,,/ - 2 (® + © re/ ) X a - 7 ~ X a?P« [ v M x] (Vg - v M ) , (29) 

where v B is the measured vector in the body and we’ve made use of the fact that the observed value of z is 
identically zero for our measurement models. Since the reference quaternion q ref is specified in the MEKF 
by the requirement that a be identically zero, d) ref is given as 

(D r ' f =m + T-'a?P a (\ B x\ M ). (30) 


The reference quaternion propagation including the effect of measurements is given by Eq. (12) with this 
angular rate. We see that the unity norm of the reference quaternion is maintained exactly by this 
propagation. No “reset” operation, as required in the discrete measurement case [18], is needed in the 
continuous measurement approximation. We emphasize that it is not necessary to integrate Eq. (29) in the 
MEKF; this equation is only used to derive Eq. (30). 

Substituting F, G , O) ref , and H M into Eq. (22) for the covariance matrix gives 

P a = ~[(ox]P a -P a [d>x] T + <t*/ 3X 3 -x-'afp( jv w | 2 / 3x3 - v„v^)p u 

( 31 ) 

- -K~ l <r; 2 {[ P a (v 8 x y M )x]P a +P a [P a (v 8 x v M ) xf }. 

The last term in this equation depends on the measurement residuals. Since the covariance is defined as the 
expectation over the distribution of measurement errors, which are assumed to have zero mean, this term 
should be omitted, giving 

P a = -[mx]P„ - PJioxf + a 2 I 3x3 -T-'a; 2 P a (|v M | 2 / 3x3 - v M v^)p a . (32) 

If this argument isn’t convincing, we can simply note that the discarded term is negligibly small compared 
to the terms retained in Eq. (32). 

The state propagation of Eqs. (12) and (30) and the covariance propagation of Eq. (32) are the results we 
need from the MEKF. 
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Unit Vector Measurements 


The errors of unit vector measurements are not isotropic, so Eq. (28) no longer applies. Consider a 
distribution of unit vectors about the z axis: 



X 

■ 

X 

u = 

y . 

= 

y 


i 

M 

1 

1 


} (x 2 + y 2 )- i(* 2 + y 2 ) 2 + ■■■ 


(33) 


We’ll assume that x and y both have Gaussian pdfs with zero mean and variance cr 2 . This is an idealization, 
since both x and y have magnitude less than unity, but it is a reasonable approximation if <x« 1. Now 


u= £{u) 


0 

0 


1-cr 2 - a* +0(a 


(34) 


and so the measurement covariance for unit vector measurements is 


= £{(u-u)(u- u) 7 } = E 


xy 


xy 

2 


<4 l-x 2 -y 2 y^ 


2 2 
’X -y 


x^l - x 1 -y 2 

yy 


vjl-x 2 - 

l-x 2 -y 


- uu 


(35) 


= a]{h+i ~ »o u o + + 0{a])] u 0 uj }, 


where u 0 = [00 if denotes the value of the unit vector for zero errors. The inverse of this is, keeping only 
the dominant term in each of the two orthogonal projections u 0 uj and / 3x3 - u 0 u£ , 


R- 


G ?V 3 


U n Un +a ' 2 


u o u o) ■ 


(36) 


When this is substituted into Eq, (21) or (22), using \ M for u 0 , only the first term in the parentheses on 
the right side of Eq. (36) gives a non-zero contribution because of the form of H M in Eq. (27). Thus unit 
vector measurements give the same result in the MEKF as vector measurements with isotropic errors. 

QUATERNION ESTIMATION 

It is natural to define the estimate of an unconstrained quaternion as an expectation value, defined as the 
integral over four-dimensional Euclidean space E4 weighted by a probability distribution function (pdf) 
p(q | Z) conditioned on the measurements Z: 

? = E{q | Z} s p(<7 | Z) d 4 <7 . (37) 


Figure 1 gives several examples of possible quaternion distributions in four-dimensional Euclidean space. 

Of course, only a two-dimensional cross-section of E4 can be shown in the figure, and one quadrant of such 
a cross-section is displayed in each plot with the origin of E4 at the lower left-hand corner. 

Unit quaternions reside on the three-dimensional sphere S3 embedded in four-dimensional Euclidean space 
E4, indicated by the arcs in Fig. 1. Thus the pdf for a distribution of unit quaternions is given by 

p(q\Z) = p(q\Z)S(\q\ 2 -l), (38) 

where p{q\Z) is the pdf on S3 and the Dirac delta function enforces the unit norm constraint. Equations 
(37) and (38) provide an unsatisfactory definition of q , however, since restricting the pdf in quaternion 
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space to the surface of a unit sphere means that its expectation value must be inside the sphere, as 
illustrated in Fig. la, unless the pdf is concentrated at a point. The unity norm constraint violation will be 
on the order of the variance of the attitude errors, like the unit vector expectation of Eq. (34). The definition 
of the covariance of a distribution restricted to S3 raises similar problems. 

It can be argued that the expectation value should not be naively performed in Cartesian coordinates in E4 as 
in Eqs. (37) and (38), but rather in terms some non-Euclidean parameterization of S3. In fact, this is exactly 
what the MEKF effectively accomplishes, although the MEKF entirely avoids the concept of a quaternion 
pdf. Having already presented this approach, we will not discuss quaternion pdfs restricted to S3 any further. 

An alternative is not to assume the quaternion to be normalized. In this case an orthogonal attitude matrix 
is formed by using the normalized quaternion q/\q\ in the quadratic form of Eq. (2), now indicated explicitly 
by a subscript Q , 

A r (q) = A Q (q/\q\) = |^|' 2 {(g 4 2 - |q| 2 )/ 3x3 + 2qq T - 2 q A [q x]j . (39) 


The subscript R identifies this as the ray representation model, since any quaternion along a ray in E4 (a 
straight line through the origin) represents the same attitude in this representation, with the exception of the 
zero quaternion at the origin. The twofold ambiguity of the unit quaternion representation corresponds to the 
two points where the ray pierces the unit sphere S3. This was referred to as the “linearized orthogonalized 
matrix” (LOM) model in References 14 and 15. Note that Eq. (5) does not hold for a non-unit quaternion. 

The 4x4 quaternion covariance matrix P q reveals its structure with the aid of the matrix 


which obeys the identities 


E(«7); 


<?4^3x3+[qx] 

-q T 


E T (q)q = 0, 


(40) 


(41) 


and 


~ T (q)~(q)=\q\ 2 I 3x3 , 

(42) 

Z(q)Z T (q) + qq T =|q| 2 / 4x4 . 

(43) 


These identities mean that the 4x4 matrix \q\ ‘ [Z(^) q\ is an orthogonal matrix whether or not q is a 
unit quaternion, so the quaternion covariance can be parameterized without loss of generality as 


p q = E ( (<? - <?)(<7 -q) T ) = [H(4 ) 


x p 

4 r a 

1J 

.2 Pc 



(44) 


The interpretation of this partitioning is that P a is the 3x3 attitude covariance, p n is the quaternion norm 
variance, and p,. is the vector of covariance between the attitude and norm errors. This factorization is not 
used in unconstrained quaternion estimation, but gives insight into what is happening behind the scenes. 

Figure lb shows a distribution for which the mean has unit norm and for which the scatter in norm (i.e. 
along the ray from the origin through the mean) is less than the scatter in angle. Figure lc shows a case for 
which the angular errors are smaller than the norm errors, and Fig. Id illustrates a pdf with isotropic errors. 
Since the expected value of the quaternion is not required to be on S3 in a ray representation. Figs, lc and 
Id illustrate cases where the norm is q = 1.1. All the pdfs illustrated in Fig. 1 have p c = 0, with the 
exception of Fig. le, which is to say that is the only figure exhibiting a correlation between norm and 
angular errors. 
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(e) (f) 

Figure 1 Representative Quaternion Probability Distributions: 

The arc is one quadrant of the projection of unit sphere S3 onto the plane. 
The centroid of the distribution in each subplot is denoted by +. 
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Finally, Fig. If shows a pdf restricted to a three-dimensional manifold in E4, similarly to Fig. la. 
However, this manifold in Fig. If is not S3, but rather is the hyperplane tangent to S3 at q . The 
covariance of the pdf illustrated by this figure has both p c = 0 and p„ =0, since all the errors are 
orthogonal to q . This also means that that q is a null eigenvector of P q . This may seem like a peculiar 
pdf, but it is the form implicitly assumed in Reference 12, and one of the EKFs we consider below will be 
shown to drive the error distribution to this limit. 


Quaternion propagation with the kinematics model of Eq. (16) can be written as 

? = + nL)®? = ^Q.((a)q + ^E(q)n m , 


with E (q) given by Eq. (40) and 


£2(<n) = 


-[cox] 

— co r 


CO 

0 


(45) 


(46) 


Equation (45) is independent of quaternion normalization, because Eq. (8) is a homogenous in q . Letting 
A q = q -q and linearizing gives the EKF kinematics 

Aq = F Aq + Gn m , (47) 

with 

F = ^Q(co) (48) 

and 

G = 2 “(<?) = 2 E (49) 

Substituting Eq. (44) into the quaternion propagation equation 

P q = FP q + P q F r + GQG t - P q H T (tRT 1 HP q . (50) 

and performing some algebra using the quaternion kinematics and Eq. (Al) of Reference 12, 


gives 


- p 

4 r a 

2 Pc 


ip! 



Pn . 

_ 


E = -^Q(ffi)E + S[ri»x], 


-±([mx]P. +Pjmx] r ) + |<r* 




t(fflxp c ) 


\P. 

2 P, 


« 2PrjH T 

I Pn If T 


~2“XPr 


\H t {zRY 


H 


1* * 


iP c 


iP. 


Pn 


(51) 


(52) 


We now consider three implementations of the quaternion EKF. 

Ray Representation EKF 

This uses the orthogonal attitude matrix of Eq. (39) in the measurement model 

z = h R (q) + n R =v B -A K (q)\ l +n z . (53) 

After a good deal of algebra, we find the measurement sensitivity matrix 

h r = 3h */ 5 <7|j =- 2 |?| 2 [v*x]S r , (54) 


with v R denoting the measurement in the body frame predicted by the expected quaternion: 

Vr = Ar( f)v,. 


(55) 
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Isotropic vector measurements with measurement covariance obeying Eq. (28) give 


E r l 

—J 1 

H t r {tR)-'H r \ e |] = 4t- , «t: 2 | 

N 2/ 3XJ-V,vj 

0 

u \ 


0 T 

0_ 


which leads to the covariance propagation equations 

K +aj / 3x3 -r-' / 3 X3 -y&)p„> 

p. = -®xp, -t‘ , <t; j p ( ,(|v j ,| j / 3x3 - vjv^Jp,, 

^ x pj 2 . 


( 56 ) 


(57a) 

(57b) 

(57c) 


The 3x3 attitude covariance P a obeys the same equation as in the MEKF, Eq. (32). We see from Eq. (44) 
that p c (f) is zero initially if the initial estimate q is an eigenvector of the initial covariance; this includes 
the case that the initial P qq is a multiple of the identity. In this case, Eqs. (57b) and (57c) show that p c (f) is 
zero for all time and the norm variance p„ is constant. The norm variance is mixed in with the attitude 
covariance in P m , which may cause loss of significance of the latter if the attitude estimate converges to 
small errors and the probability distribution tends to be concentrated along a ray, as in Figure lc. 


The expectation of Eq. (45) augmented by a measurement term gives the equation for the time dependence of 
the quaternion estimate 


f = ±Q(m)f-PffJ(T RT'[ v fl 


- Ar (<?) v, ] = iQ( 0 >) q + 2t "‘a; 2 


E q\ 


x p 

4 r a 

iPl 


ip, 

Pn 


[v«xf 

o r 


(58) 


where we have made use of the fact that [v s xf A R {q)v, ~ -v R xv,=0. This can be expressed as 

f = ® f + (▼* X p e ) ' V fi ] f. 


with 


''^.(v.xv,). 


(59) 

(60) 


The norm of the quaternion estimate satisfies 


d_ 

dt 



(61) 


We see that the norm is constant if \ R x (r) = 0, which is consistent with Eq. (57c). The time 
dependence of the normalized quaternion is found by combining Eqs. (59) and (61): 


f 

d_ 

dt 

V 



q_ 

|?| 


q d 




4 f 


®q . 


( 62 ) 


Thus the attitude estimate is the same for the ray representation EKF as for the MEKF. 

Using Eq. (36) for unit vector measurement errors gives, remembering that v R = u 0 is a unit vector, 



which is the same as the form derived above for vector measurements with isotropic errors. Thus the ray 
representation EKF gives the same covariance and state updates with unit vectors as with vectors having 
isotropic errors. 

The ray representation form of quaternion EKF, which we have seen to be effectively equivalent to the 
MEKF, has been applied to attitude estimation of the ALEXIS and CAPER spacecraft [29, 30]. 

Quadratic Representation EKF 

This uses the homogenous quadratic expression for the attitude matrix of Eq. (2) in the measurement model 

z = h Q (q) + n Q =y B -A Q (q)\, + n.. (64) 

We note that the attitude matrix is not orthogonal unless the quaternion has unit norm. We will see that 
these measurements drive the norm to unity, however. The measurement sensitivity matrix is 

H q =dh Q fdq\^ = —2 !^| ([v e x]E T + v Q q T ) , (65) 


where 

v Q = A Q (q)\, 

is the predicted measurement in the body frame. For isotropic measurement errors, we have 


1 

—T j 

H T Q (zRy'H Q [s f ] = 4t~ i <7' 2 


0 

i i 2 

L ? J 

i 

o T 

M . 


which leads to the covariance propagation equations 


P a = -[tox]P a - PJwxf + crX 3 -*"<T 


-1 _ -2 


n n /m-vS j /> a + 4|v e i Pc p ' 


P r =4x Pc -r-'cr; 2 


'Q’e 

P a\ \v Q \ / 3 *,-v c v' j + 4|v e | p n 

|2 


Pc 


Pn =-^Xl |^<2 X Pc| +4 K| Pi ■ 


( 66 ) 

(67) 

(68a) 

(68b) 

(68c) 


The last term in Eq. (68c) drives p„ asymptotically to zero, and the last term in Eq. (68b) drives p c to zero. 
Thus the 4x4 covariance matrix becomes singular asymptotically with null vector q , and the probability 
distribution approaches the form illustrated in Figure If. The resultant potential instability of the quadratic 
representation EKF can only be avoided by adding unphysical process noise to the quaternion magnitude. 

The propagation equation for the state estimate is 

% = ±Q(®)f + PH t q (r R)-'[ z Q - h 2 m = i£2(m) f - PH T Q {rR Q T\y B - A e (f) v, ] . (69) 


This gives 


which can be written as 


with 


f- -i 

i P a 2P f 

"[v e xf' 


>P c Pn . 

v r 

L v e J 


(v B -v e ), 


<2 = 2 ® <7 +T‘'a; 2 [(v e xp c + 2 p„y Q ) ■ (v B - y Q )]q, 


s to + z-'a; 2 [P a (v s x v c ) + 2[v e ■ (v a - v fl ) ]p e } . 


(70) 

(71) 

(72) 
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The quaternion norm satisfies 


d_ 

dt 


= 2 q T q = 2r' l OV 


i 2 

<7 


(v c xp f . + 2/?„v e )-(v s 



(73) 


If q is not a unit quaternion, \ B - v Q will have an error along \ Q that gives rise to a correction driving 
the norm of q toward unity as long as p n is not zero. If p ( . (r) is zero initially (and thus for all time), the 
angular behavior given by Eq. (72) is the same as for the ray representation EKF and MEKF; but Eq. (71) 
shows a non-zero norm update parallel to q , as is expected from the norm variance. If p c (t) is not zero, 
the attitude behavior will also differ from the ray representation EKF and MEKF, and it is difficult to avoid 
the conclusion that the quadratic representation EKF will give incorrect results in this case. 

For unit vector measurements, the predicted body vector \ Q is not a guaranteed to be a unit vector in the 
quadratic measurement EKF. However, \ Q is parallel to u 0 , the unit vector in Eq. (36), so 


H C R ml< = ~ 2a T k (H [ V G X f + or ; 2 q ) . 


(74) 


This differs from the result for vector measurements with isotropic errors 


HqR~ 1 =-2 a: 


k (“[VnXf +qy T 0 ) ■ 


(75) 


The result is to multiply the (4,4) element of the matrix on the right side of Eq. (67) by a~ 2 for unit vector 
measurements. Both p n and d\q | /dt gain an extra factor of cr~ 2 when p, is zero. Since a z can be very 

small, the approach to quaternion normalization and also to covariance matrix singularity can be much more 
rapid for unit vector measurements than for measurements of vectors with isotropic errors. The explicit time 
dependence of p„ for p c (t) = 0 and jv^] constant is 


p„('H 


p;'(0) + 4a: 2 |v c | 2 (t/T) 
p; 1 (0) + 4a; 4 |vJ 2 (t/T)] 1 


for vectors with isotropic errors 
for unit vectors. 


(76) 


Note that the convergence of |<y| to unity and of p n to zero is a result purely of measurement processing. 

Neither “brute force” normalization of the quaternion nor norm-enforcing pseudo-measurements is required, 
in principle. Such subterfuges have been found useful in the presence of norm errors arising from numerical 
and discretization effects in actual implementations, however [13-15]. 


The Linear Measurement EKF 


Choukroun, Bar-Ithack, and Oshman formulated a measurement model that is linear in the quaternion [31]: 

Z = h z.(^) + n z =2(^)v fi -¥(;?) v, +E(^)n ? , (77) 

where 





(78) 


Equations (41)-(43) and the identity 

E T (q)V(q) = A Q (p) = \q\ 2 A R (q) (79) 


indicate that this model is closely related to the ray representation EKF. 
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It is easily seen that 


h L (q) = Z(q)\ B - x ¥(q)\, 


-[(Vb+vJx] (Vg-V,) 


-( v a-v ; ) 


0 


q=H L q. 


(80) 


In practice, we must use the predicted measurements to compute the sensitivity matrix. Substituting 
v fl = v R and v, = Al(q)v R into the above equation gives, after some algebra. 


'E[v R x]E r =Etf*. 


(81) 


We also use q in the measurement noise covariance, giving the singular measurement covariance 

R l =E£{nV)E r = cr 2 SS r =cr;lim(EE r +eff r ). (82) 

s —*0 

However, Eqs. (41)— (43) and (81) give 


IT T r>-l _ tjT ~T 
** L ~ "R i ~ 


CT; 2 |<5f|' 4 lim(EE r +e 


■' qq T ) 


= a 


ir^s r - 


(83) 


Comparing with Eqs. (28) and (54) shows that this model gives the same covariance and state propagation 
as the ray representation EKF. It is not difficult to see that this equivalence also holds for unit vector 
measurements. The apparent simplicity of the linear dependence of the measurements on the state vector 
given by Eq. (77) is compensated by the appearance of state-dependent measurement noise, and the result is 
mathematically equivalent to ray representation EKF with its nonlinear measurement model. 


DISCUSSION 


Of the four methods considered in this paper, the MEKF requires the least computational effort due to the 
lower dimensionality of its state vector and covariance matrix. The MEKF is also the most satisfying 
conceptually, since it respects the dimensionality of the rotation group, and its attitude estimate is a unit 
quaternion by definition. 

The ray representation EKF and the equivalent linear measurement EKF give the same attitude covariance 
and attitude estimates as the MEKF, at the cost of some additional computational burden. They are 
conceptually simpler than the MEKF, since they require no attitude parameterization other than the 
quaternion. Some numerical significance issues may arise from the unobservable norm degree of freedom of 
the quaternion, but these have not caused difficulties in practical applications. 

The quadratic representation EKF rests on the least secure foundation. We have shown that the quaternion 
approaches unit norm, and thus the attitude matrix approaches an orthogonal matrix, as a natural result of 
measurement processing, without resorting to “brute force” normalization of the quaternion or norm- 
enforcing pseudo-measurements. However, the quaternion norm variance and attitude-norm correlations tend 
asymptotically to zero for this filter, so the 4x4 covariance matrix becomes singular. This singularity, 
which may lead to instability of the Kalman filter, can only be avoided by adding unphysical process noise 
to the quaternion magnitude. The singularity is approached more rapidly for unit vector measurements than 
for vector measurements with isotropic errors. 

If the initial quaternion estimate is an eigenvector of the initial 4x4 quaternion covariance in the quadratic 
representation EKF, the ray representation EKF, or the linear measurement EKF, the correlation between 
the attitude errors and quaternion norm errors will be zero initially and for all time. For this initialization, 
the time dependence of the attitude covariance in these three filters is the same as the MEKF, independent of 
the quaternion norm variance. If the correlation between the attitude errors and quaternion norm errors is not 
zero, the attitude covariance of the quadratic representation EKF differs from that of the other filters. Also, 
the attitude matrix provided by quadratic representation EKF is not guaranteed to be orthogonal, so the 
estimates provided by this method should be regarded with caution. 
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The motivation for considering the four-component quaternion EKFs despite the conceptual and 
computational advantages of the MEKF appears to be the superficial resemblance of a quaternion EKF to a 
linear Kalman filter. This resemblance is deceiving, because the process noise and any dynamic parameters 
to be estimated enter mutiplicatively rather than additively in the quaternion kinematics equation. Most 
measurement models are also inherently nonlinear, since we have seen that the apparent simplicity of the 
linear measurement model is compensated by the appearance of state-dependent measurement noise, and the 
result is mathematically equivalent to ray representation EKF with its nonlinear measurement model. Thus 
we see no valid reason to prefer a quaternion estimator to the MEKF. 
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